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Abstract 

The basic problem in equilibrium statistical mechanics is to compute 
phase space average, in which Monte Carlo method plays a very important 
role. We begin with a review of nonlocal algorithms for Markov chain 
Monte Carlo simulation in statistical physics. We discuss their advantages, 
applications, and some challenge problems which are still awaiting for 
better solutions. We discuss some of the recent development in simulation 
where reweighting is used, such as histogram methods and multicanonical 
method. We then discuss the transition matrix Monte Carlo method and 
associated algorithms. 

The transition matrix method offers an efficient way to compute the 
density of states. Thus entropy and free energy, as well as the usual 
thermodynamic averages, are obtained as functions of model parameter 
(e.g. temperature) in a single run. 

New sampling algorithms, such as the flat histogram algorithm and 
equal-hit algorithm, offer sampling techniques which generate uniform 
probability distribution for some chosen macroscopic variable. 

1 Statistical Physics 

Statistical mechanics is a theory developed at the end of nineteenth century 
to deal with physical systems from an atomistic point of view. In principle the 
properties of bulk matter, which may contain 10^^ atoms, can be worked out 
from the motion of atoms following the basic equations of Newtonian mechanics 
or quantum mechanics. However, such detailed information is not available or 
not really necessary. A probabilistic point of view, in the form of statistical 
ensemble is taken. The theory is extremely economical and successful in dealing 
with equilibrium problems. 

There are a number of equivalent formulations of the theory. But statistical 
mechanics in a nutshell is the following concise formula that connects statistical 



mechanics with thermodynamics. First the partition function is defined as 

Z = ^ exp 

X 

where X is a "state" of the system; the summation is carried over ah possible 
states. T is absolute temperature and k is the Boltzmann constant (1.38 x 10"^"^ 
Joule/Kelvin). The free energy of the system at temperature T is given by 

F = -kT\TLZ. (2) 

Free energy is a useful thermodynamic quantity in dealing with phase transi- 
tions. Other macroscopic observables are averages of the corresponding micro- 
scopic quantities over the Boltzmann-Gibbs weight exp(— i?(X)/fcr), 

{9{X)) = ^Y.^{X)e^^{-E{X)/{kT)). (3) 

X 

The following three quantities are perhaps the most important ones, internal 
energy U ^ (E), heat capacity C ^ dU/dT = {{E"^) - {E)'^) /{kT'^), and entropy 
S = {U — F)/T = —k{\np{X)) , where the average is over the Gibbs probability 
density, p{X) = c^p{-E{X)/kT)/Z. 



2 Monte Carlo Method 

The essential task in statistical mechanics is to do the multi-dimensional in- 
tegrations (for continuous systems) or summations (for discrete systems). To 
fix the notation and language, we briefly introduce the basics of Monte Carlo 
method ^, ^. Consider the computation of a statistical average, 

G= f g{X)p{X)dX={g), (4) 

where the probability density obeys p{X) > 0, and J^p{X)dX = 1. Suppose 
that we can generate samples Xi according to the probability p{X), then the 
integral can be estimated from an arithmetic mean over the samples, 

1 

GM^^Y.9iX^)■ (5) 



The random variable Gm has mean G and standard deviation y/T{{g — G)^) /M, 
T being decorrelation time. Thus in the limit of large number of samples, the 
estimate converges to the exact value. 

The most general method to generate X according to p{X) is given by a 
Markov chain [|j with a transition probability P{X'\X) = W{X X'), satis- 
fying the conditions stated below. This is the probability of generating a new 
state X' given that the current state is X. Such process converges if 

W"{X' X) > 0, for aU X' and X, and n > no. (6) 
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The equilibrium distribution of the Markov chain satisfies 

= (7) 

X' 

In constructing a Monte Carlo algorithm, it is convenient to consider a much 
stronger condition, the detailed balance 

p{X')W{X' ^ X) =p{X)W{X ^ X'). (8) 

One of the most famous and widely used Monte Carlo algorithms is the 
Metropolis importance sampling algorithm Q ^ . It takes a simple choice of the 
transition matrix: 

W{X ^ X') ^ Six ^ X') min (^1, ^) (9) 

where X ^ X', and is a conditional probability of choosing X' given that the 
current value is X, and it is symmetric, S{X X') — S{X' X). Usually 
S{X X') = unless X' is in some "neighborhood" of X. The diagonal term 
of W is fixed by the normalization condition J2x' ~^ — 1- 



3 Cluster Algorithms 

In order to facilitate the discussion we first introduce the Ising model. Ising 
model is an interacting many-particle model for magnets. A state consists of a 
collection of variables ai taken on only two possible values +1 and —1, signifying 
the spin up and spin down states. The spins are on a lattice. The energy of the 
state (T = {ci} is given by 

S((t) = -J^CT.CTj-, (10) 

where J is a constant which fixes the energy scale, the summation is over the 
nearest neighbor pairs. When temperature T is specified, the states are dis- 
tributed according to 

p{a) = Z-^e.p^-^y (11) 

In a local Monte Carlo dynamics (Metropolis algorithm), one picks a site i 
at random, i.e., choosing a site with equal probability (this specifies or realizes 
S). Then, the energy increment if the spin is flipped, <Ti — > — ct;, is computed, 
which has the result AE = 2JaiJ2j (^j- The flip is accepted with probability 
min(l, exp(— Ai?/fcT)) . If the flip is rejected, the move is also counted and the 
state remains unchanged. One Monte Carlo step is defined as N moves (trials) 
for a system of N spins, such that each spin is attempted to flip once on average. 
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The local algorithm of Metropolis type has some salient features: (1) it is 
extremely general. Less assumption is made of the specific form of probability 
distribution. (2) Each move involves 0(1) operations and 0(1) degrees of free- 
dom. (3) The dynamics suffers from critical slowing down. The correlation time 
T diverges as a critical temperature is approached. We shall elaborate more on 
this in the following. 

The statistical error, using estimator Eq. (a), is given by 



where Var(G'i) = cr^ = (g^) — (5)^ is the variance of the observable, M is the 
number of Monte Carlo steps, and t is the decorrelation time. We can take the 
point of view that the above equation defines r, i.e., 

MVar(GM) 
Af^oo Var(G'i) ^ ' 

Perhaps it is appropriate to call r decorrelation time, since (r — l)/2 is some- 
times called correlation time in the literature. The decorrelation time r is the 
minimum number of Monte Carlo steps needed to generate effectively indepen- 
dent, identically distributed samples in the Markov chain. The smallest possible 
value for r is 1, which represents independent sample by every step. The usual 
integrated autocorrelation time j|] differs from our definition by a factor of 2, 

Tint — T /2. 

The critical slowing down manifests itself by the fact that t (x at the 
critical temperature Tc where a second order phase transition occurs. Here 
L is the linear dimension of the system {N — L'^ in d dimensions). For the 
local algorithms for many models and in any dimensions, z « 2. This suggests 
bad convergence, specially for large systems. At a first-order phase transition, 
where some thermodynamic variables change discontinuously, the situation is 
even worse — r diverges exponentially with system sizes. 

For the two-dimensional Ising model, a phase transition occurs |^ at kTc/ J « 
2.269. The magnetization is non-zero below this temperature and becomes zero 
above Tc- In addition, in the limit of large system (L — > 00), heat capacity per 
spin and fluctuation of magnetization diverge. These intrinsic properties make 
computer simulation near critical point very difficult. 

Cluster algorithms |l0| overcome this difficulty successfully. For example, 
for the two-dimensional Ising model, the dynamical critical exponent defined in 
r oc is reduced from 2.17 for the single-spin flip to much small value for 
the cluster algorithm |l^, |l^, |l^ . It turns out that a precise characterization of 
the Swendsen-Wang dynamical critical exponent is very difficult, due to weak 
size dependence of the decorrelation time. Table 1 represents a recent extensive 
calculation, based on the definition Eq. ( ]l3| ) for the total energy, rather than the 
usual method of extracting information from time-dependent auto-correlation 
functions. In this calculation, the variance of the sum of M consecutive energies 
(c.f. Eq. (|^)) arc computed explicitly. Good convergence to the limiting value 
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L 



T 



standard error 
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4.04575 

5.17696 

6.5165 

8.0610 

9.794 

11.723 

13.872 

16.29 

18.87 



0.00033 

0.00032 

0.0012 

0.0018 

0.004 

0.012 

0.009 

0.04 

0.2 



8 

16 
32 
64 
128 
256 
512 
1024 



Table 1: The decorrelation time r of the Swendsen-Wang dynamics of the two- 
dimensional Isng model for different linear lattice size L at the critical temper- 
ature. 

is already achieved typically for M w 10^. An \/M extrapolation is used to 
get more accurate estimates of the limit. From this calculation, we find for 
the two-dimensional Ising model, the convergence to asymptotics are slow. It 
appears that the divergence is slightly faster than logL. If we fit to power law 
L^"™, using two consecutive sizes of L and 2L, the exponent Zsw decreases from 
0.35 to 0.21. It is not clear whether Zsw will convergence to a finite value or 
continue decreasing to in the limit of i — > cxd. 

The nonlocal cluster algorithm introduced by Swendsen and Wang for the 
the Ising model (and more generally the Potts model) goes as follows: (1) go 
over each nearest neighbor pair and create a bond with probability 



That is, if the two nearest neighbor spins are the same, a bond is created between 
them with probability 1 — exp(— 2 J/ZcT); if spin values are different, there will 
be no bond. (2) Identify clusters as a set of sites connected by zero or more 
bonds (i.e., connected component of a graph). Relabel each cluster with a fresh 
new value -1-1 or —1 at random. 

We note that each Monte Carlo step per spin still takes 0(1) in computa- 
tional cost. The method is applicable to models containing Ising symmetry, i.e., 
the energy is the same when at is changed to —cji globally. 

The algorithm is based on a mapping of Ising model to a random cluster 
model of percolation. Specifically, we have |]l^ |l^ |l^, |l^ 




(14) 



Z 




(15) 



{^} ('-j) 

XI XI n [p^'^.^'^AmjA + (1 - P¥n,,A 



(16) 



{^} {"} (^i) 




(17) 



{"} 
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where K — 2J/kT, p = 1 — exp(— if), Sij is the Kronecker delta, and b is 
the number of bonds, d is the dimension of a simple hypercubic lattice, and 
A^c is the number of clusters, ct; is spin on the site and riij = 0, 1 is the bond 
variable between the sites i and j. It is evident that the moves in the Swendsen- 
Wang algorithm preserves configuration probability of the augmented model 
containing both the spins and bonds. 

A single cluster variant due to Wolff Q is very easy to program. The 
following C code generates and flip one cluster: 

#define Z 4 
#define L 8 
double p = 0.7; 
int s [L*L] ; 

void flip (int i, int sO) 
{ 

int j , nn [Z] ; 

s[i] = - sO; 
neighborCi, nn) ; 
for(j =0; j < Z; ++j) 

if(sO == s[nn[j]] && drand48() < p) 
f lip(nn[j] , sO) ; 

> 

In this single cluster version, a site i is selected at random. The value of the spin 
before flip is sq. It flips the spin and looks for its neighbors. If the value of the 
neighbor spins are the same as so, with probability p a neighbor site becomes 
part of the cluster. This is performed recursively. It turns out that single cluster 
is somewhat more efficient than the original Swendsen-Wang, particularly in 
high dimensions. For dimensions greater than or equal to 4, Swendsen-Wang 
dynamics gives the dynamic critical exponent z = I, while Wolff single cluster 
is z = 

It is easy to see why the single cluster algorithm of the above works. Let 
us consider a general cluster flip algorithm with two bond probabilities: Ps will 
be the probability of connecting two parallel spin sites; Pd the probability of 
connecting anti-parallel sites. Consider the transition between two configura- 
tions A and B, characterized by flipping a cluster C. A cluster is growing from 
a site i until the perimeter of the cluster is not connected to the outside. The 
transition probabilities can be written down as 

wiA^B)^pj n (1-^^) n (1-^^)' (18) 

ac(TT) dciu) 

w{B^A)^pj n i^-Pd) n (^-p^y (19) 

ac(TT) dciu) 
Note that the bond configuration probabilities are the same in the interior of 
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the clusters. The difference occurs at the boundary dC, where parallel spins 
(tt) in ^ becomes anti-parallel spins (ti) in B, or vice versa. Detailed balance 
requires that 

WiA -^B) _ P{B) „ 



W{B A) P{A) 



(1 - P.y^-'^'il - p^)^v2-ivi ^ g— ^ (20) 



where A'^i is the number of parallel spins on boundary of C in configuration A; 
and A^2 is the number of anti-parallel spins on boundary of C in configuration 
A. Since we have AE = Eb - Ea = 2{Ni - N2)J, we obtain 

Although the algorithm is valid for any < < 1, it is most efficient at P^; = 0, 
the Coniglio-Klein bond probability value. 

A quite large number of statistical models can be treated with cluster algo- 
rithms, with varied success. Excellent performance has been obtained for Ising 
model, Potts models, and antiferromagnetic Potts models XY model and 
general 0{n) models 0^ field-theoretic model |2|], some regularly frustrated 
models ||2^, , six- vertex model [|2^ , etc. Cluster algorithms are proposed for 
hard sphere fluid systems [ p5[ , quantum systems p6| , [2^ , p8| , p9[ , microcanon- 
ical ensembles conserved order parameters ||3l|, |32|, etc. Invaded cluster 



algorithm 1 33 and other proposal |Q are excellent methods for locating critical 
points. The cluster algorithms are also used in image processing |3^ . 

The cluster algorithms do not help much in temperature-driven first-order 
phase transition - the slow convergence has been shown rigorously p^ . Mod- 
els with frustration, spin glass being the archetype, do not have efficient cluster 
algorithms, although there are attempts |38[|9l |§ with limited success. Break- 
through in this area will have a major impact on the simulation methods. 



4 Reweighting Methods 

In this and the following sections, we discuss a class of Monte Carlo simulation 
approaches that aim at an efficient use of data collected, and sampling methods 
that enhance rare events. 

The computation of free energy, Eq. (||) , poses a difficult problem for Monte 
Carlo method. A traditional method is to use thermodynamic integration, e.g., 

P(T.) FiT,)_ r^iEl^^^ ^^^^ 



T2 Ti Jt, T2 



based on the relation (E) = —d\nZ/df3, where (3 — l/(kT). 

If we can estimate the density of states (the number of states with a given 
energy E for discrete energy models), then we can compute free energy, as well as 
thermodynamic averages. The result is obtained as a function of temperature 
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T, rather than a single datum point for a specific value of T, as in standard 
Monte Carlo simulation. 

This idea has been pursued over the last decade by Ferrenberg and Swendsen 
H, Berg et al g, |4|], Lee |6|, Oliveira et al @, and Wang |§. 
Consider the following decomposition of summations over the states 



^A(a)ei-^) .= ^6-^ ^ A(a)=5]e-^n(i?)(A(a))^, (23) 

{<t} E E{a)=E E 

where {A) e is the microcanonical ensemble average, 

^ ' E{cr)=E 

Since the state space is exponentially large (2^ for the Ising model with N 
spins), and the range of E is typically of order N, if n{E) can be computed 
accurately, the task is done. The canonical average of A is related to the micro- 
canonical average through 



EE{A)Ee^v{~^)n{E ) 
E^exp(-J.)n(£;) 



iA)r = ^^-/^;;-^v^.t;-v-^ ^ (25) 



and free energy is computed as 



F = -kT\nJ2e^p(^~^^n{E). (26) 
4.1 Histogram method 

Ferrenberg and Swendsen popularized a method which in a sense is to 
compute the density of states (up to a multiplicative constant) in a range close 
to a given simulation temperature. This method is generalized as multiple 
histogram method to combine simulations at differential temperatures, to get 
the whole energy range |Q . We discuss here only the single histogram method 
for its simplicity. 

During a normal canonical simulation at fixed temperature T*, we collect 
the histogram of energy, H{E), which is proportional to probability distribution 
of energy, 

H{E)=cn{E)cyiY>{-E/{kT*)). (27) 

The constant c is related to the partition function, c — M/Z{T*), where M is 
the total number of samples collected. From the above equation, we find n{E) cx 
H{E)exp{E/{kT*)). With this information, we can compute the free energy 
difference between temperature T* and a nearby temperature T. Similarly, 
moments of energy can be computed after the simulation, through histogram 
reweighting, 

_ j:EE"H{E)c^p{-E/ikT) + E/ikT*)) 

EEH{E)c^v{-E/{kT) + E/{kT*)) ' ^ ^ 
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The range of E that the histogram data can be collected at a fixed temper- 
ature is limited by the energy distribution, which for the canonical distribution 
away from critical point, is of order of ^/N. The whole range of energy E is of 
order N. This limit the usefulness of single histogram method. 

4.2 Multicanonical Monte Carlo 

The multicanonical Monte Carlo method has been shown to be very effective 
to overcome supercritical slowing down, reducing the relaxation time from ex- 
ponential divergence with respect to system size to a power, at the first-order 
phase transitions [Q. Multicanonical ensemble flattens out the energy distri- 
bution, so that the computation of the density of states n{E) can be done for all 
values of E. A multicanonical ensemble is defined to following the probability 
density for the states as 

= f{Ei.)) oc (29) 

such that the energy histogram H{E) (cx n{E)f{E)) is a constant. From the 
histogram samples obtained by a simulation with the weight of state at energy E 
as f{E), the density of state can be computed |46| from n{E) — H{E)/ (^f{E)c). 
However, unlike canonical simulation where f{E) = e"^/'^^ is given, in a mul- 
ticanonical simulation, f{E) is unknown to start with. 

Berg proposed an iterative method to compute the weight in a parametrized 
form f{E) = exp[—f3{E)E + a{E)), starting with no information, fo{E) = 
const. A new estimate at iteration n is then based on the results of all previous 
iterations. We refer to references Q, ^ for details. 

5 Transition Matrix Monte Carlo and Flat His- 
togram Algorithm 

The flat histogram algorithm offers an efficient bootstrap to realize the multi- 
canonical ensemble, while transition matrix Monte Carlo utilizes more data that 
can be collected in a simulation to improve statistics. 

5.1 Transition matrix 

We start from the detailed balance equation for some given dynamics: 

p{a)W{a a') = p{a')W{a' a). (30) 

By summation over the states a of fixed energy E, and a' of fixed energy E' , 
and assuming that the probability of the state is a function of energy only, 
p(cr) cx f{E{a)), we get 

n{E)f{E)T{E ^ E') = n{E')f{E')T{E' ^ E), (31) 
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where the transition matrix in the space of energy is defined as 

^(^-^^'^^^ ^ ^ W^(a-a'). (32) 

^ ' E{cr)=E E{a')=E' 

The matrix T has a number of interesting properties: it is a stochastic matrix in 
the sense of T{E ^ E') > and J2e' T{E E') = 1; the stationary solution 
of T is the energy distribution n{E)f{E)\ the dynamics associated with T is 
considerably faster than that of W . 

We specialize to the case of single-spin-flip dynamics for the Ising model. 
The transition matrix W for the spin states consists of a product of two factors, 
the probability of choosing a spin to flip S{a a'), and the flip rate a{E — > 
E') = min(l, f{E')/f{E)). We have S{a -> cr') = less the two configurations 
(7 and cr' differ by one spin, in this case, the value of S is where N is the 

number of spins in the system. Using these results, we can rewrite the transition 
matrix as 

T{E^E') = ^{N{a,E' -E))^a{E^E'), E ^ E' . (33) 

The diagonal elements are determined by normalization. Substituting Eq. 
into Eq.(^l|), using the relation between a{E E') and ,f{E), we obtain 

n{E){Nia,E' - E)) ^ = n{E'){N{a\ E - E'))^,. (34) 

This is known as broad histogram equation ^ which forms the basis for the 
flat histogram algorithm presented below. Additionally, this equation also gives 
us a way of computing the density of states n{E) by the quantity {N{a, E'~E)) e 
obtained from spin configurations generated from any distribution /(£'((t)). 
The quantity N{(t, AE) is the number of ways that the system goes to a state 
with energy E + AE, by a single spin flip from state a. The angular brackets 
indicate a microcanonical average: 

{Nia,AE))^^-^ Nia,AE). (35) 

^ ' E{a)=E 

5.2 Flat histogram algorithm 

The following algorithm ^ generates flat histogram in energy and realizes 
the multicanonical ensemble. 

1. Pick a site at random. 

2. Flip the spin with probability 

,P . A {N{a',~AE))E+AE \ 

aiE^E )= mm [l, ^^(^^^^)^^ ) , (36) 

where the current state a has energy E, the new state a' has energy 
E' = E + AE. 
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3. Accumulate statistics for {N{(t,AE))e- 



4. Go to 1. 

We note by virtue of Eq. (34), the flip rate is the same as that in multicanonical 
simulation with a weight l/n{E) and Metropohs acceptance rate. While in 
multicanonical sampling, the weight is obtained through several simulations 
iteratively, the quantities {N{a, AE)) e is much easier to obtain, through a single 
simulation. This quantity serves a dual purpose — it is used to construct a Monte 
Carlo algorithm (used as input), and at the same time, it is used to compute 
the density of states (output of the simulation) . Clearly, this is circular unless 
approximation is made. We have considered replacing the exact microcanonical 
average by an accumulative average, over the history of simulation generated so 
far, i.e., 

1 ^ 

(Nia, AE))^ « J2 SEia^),ENia\AE), (37) 

where {a\ i = 1, 2, • • •} is the sequence of states generated with the algorithm 
given above; H{E) is the number of samples accumulated at the energy bin E. 
In case the data for computing the flip rate is not available, we simply accept 
the move to explore the new state. A more rigorous way of doing simulation 
is to iterate the simulation with constant flip rate. For example, after the 
first simulation, we compute a first estimate to the density of states. In a 
second simulation, we perform multicanonical simulation a la Lee The 
data collected in the second run for {N{(t, AE)) e will be unbiased. It is found 
that even with a single simulation, the results converge to the exact values 
(for {N{a, AE)) E and n{E)) for sufficiently long runs, even though a rigorous 
mathematical proof of the convergence is lacking. 

Wang and Landau js^ proposed recently a new algorithm that works directly 
with the density of states n{E). The simulation proceeds with the flip rate 
min(l, n{E) /n{E')) , but the value of the density of states is updated after every 
move by n{E) ^ n{E)f and letting / 1 for convergence. Excellent results 
were obtained. A careful comparison with flat histogram method is needed. 



5.3 N-fold way (rejection-free moves) 

In Metropolis algorithm, moves are sometimes rejected. This rejection is impor- 
tant for realizing the correct stationary distribution. In 1975 Bortz, Kalos, and 
Lebowitz [|^ proposed a rejection-free algorithm. It is still based on Metropohs 
flip rate, but the waiting due to rejection is taking into account by consider- 
ing all possible moves. The Bortz-Kalos-Lebowitz N-fold way algorithm for the 
Ising model goes as follows: 

1. Compute the acceptance probability for one attempt of a move 

^-E^^^«(^-^ + Ai.). (38) 

AE 
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2. Pick an energy change Ai? according to probability 

PAE^^^^^^a{E-^E + AE). (39) 

flip a site belonging to AE with probability 1. The site is choosing from 
the N{a, AE) sites with equal probability. 

3. One N-fold-way move is equivalent to 1/A moves in the original dynam- 
ics. Thus thermodynamic averages are weighted by 1/A, i.e., (g) = 
'^f{gt/A) I 1/^, where summation is over every move t. 

In order to implement step 2 efficiently, additional data structure is needed so 
that picking a spin in a given class characterized by AE is done in 0(1) in 
computer time. 

5.4 Equal-hit algorithm 

Combining N-fold way and flat histogram algorithm is easy, since the important 
quantity N{a, AE) is already computed in flat histogram algorithm. The flip 
rate a is given by formula (|3^). In the flat histogram algorithm, the probability 
that the energy of the system is E' is a constant, i.e., 

H{E) = {SEia),E) <^n{E)f{E)= const 
1 \ / / 1 



The averages in the second line of the above equation refer to samples generated 
in an N-fold way simulation. 

In equal-hit algorithm (ensemble) ||5^ , we require that the number of "fresh 
configurations" generated at each energy is a constant. More precisely equal-hit 
ensemble is defined by 

u{E) = {Se{^),e)n = const. (41) 
One possible choice of the flip rate is 

(P J,', ■ (, {^)E',N{Nia',E-E'))E' \ 

a{E^E) = mmll,—^ ^^T^ ' (^^) 

V {^)E.N{N{a,E' - E))e J 



where 



1 \ _ {^E{a).EA^)N 



A/ E,N {SE(a),E)N 



(43) 



is the inverse total acceptance rate 1/A arithmetic averaged over the N-fold way 
samples at energy E. 

The histogram H{E) generated in the equal- hit algorithm depends on the 
precise dynamics (the rate a) used. Since there are many possible choices of the 
rate, such "equal-hit ensemble" is not unique. 
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5.5 Determination of the density of states 

While Eq. (^J) gives us a way of obtaining the density of states, there are 
more equations than unknowns. We consider two optimization methods. The 
first method is based on the transition matrix itself. We define Te^ae = 
{N {a, AE)) E / N . Symbols with hat being unknown, and Te,ae the Monte 
Carlo estimate, consider 

Minimize ^ ^ e^,ae {Teae-Teae) , (44) 

E,AE 

subject to < Teas < 1, J2ae^eae = 1, a,iidfE.ifE+isTE+2,-2 = fE,2TE+2-iTE+i 
The last constraint needs more explanation. We assume that the energy level is 
equally spaced (as in the Ising model). Consider three energy levels, E, E + 1, 



E + 2. If we write down three equations of type (34), for transitions from E to 
E + 1, E + 1 to E + 2, and E + 2 to E, we can cancel the density of states by 
multiplying the three equations together. This leaves the last equation above, 
and it is known as TTT identity. It can be shown that multiple T identities 
(four or more) are not independent, and they need not put in as constraints. For 
Ising model there is also one additional symmetry constraint, j = T^e.-i- 

When the solution for T is found, we can use any of the energy detailed 
balance equation to find density of states n{E). The TTT identity guarantees 
that the answer is unique whichever detailed balance equation is used. 

The second method is based on optimization directly with variable n{E), 
actually S{E) = lnn{E), by 



minimize ^ a'^^ (siE + AE) - SiE) - In ^^'^^ ) , (45) 
\ J-e+ae,-aeJ 



subject to, for the Ising model, J^e "'(-^) — 2^, where N is the total number of 
spins in the system. In addition, we can put in the known fact that the ground 
states are doubly degenerate, n{0) — 2. 

Figure 1 shows one of the simulation results of the density of states, using 
the second method. The errors by comparison with exact known values |56| are 
presented in the insert of the figure. The density of states n{E) is determined 
to an accuracy of better than 2 percents in a matter of few minutes of computer 
time. 

The flat-histogram dynamics is used to study spin glasses |]5^ . The dynamic 
characteristics is quite similar to multi-canonical method of Berg. The study 
of lattice polymer and protein folding is under way. For related ideas and 
approaches, see refs.[Q ^ |60| . 

6 Conclusion 

The phenomenon of critical slowing down can be effectively dealt with by cluster 
algorithms for a large class of statistical mechanics models. We reported new 
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Figure 1: Logarithmic density of states, S{E) — lnn{E), calculated by transi- 
tion matrix Monte Carlo method for a 16 x 16 two-dimensional Ising model with 
a total 1.1 X 10^ Monte Carlo steps and the first 10^ discarded, using the flat- 
histogram algorithm and N-fold way. The insert is error in S{E) with respect 
to exact results. 



and very accurate results for the decorrelation times of the Swendsen-Wang dy- 
namics. Large size asymptotic behavior is analyzed. For super-critical slowing 
down occurring in first-order phase transitions, multicanonical ensemble simu- 
lation and flat-histogram or equal-hit algorithms arc very effective. Since the 
latter algorithms and associated transition matrix method is efficient in com- 
puting the density of states, this method can also be useful for general counting 
problems by Monte Carlo method. 
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